% Generate the data and summary statistics presented in Tables B1 and B2.

clear
clc

% Designate working directories
dir = '.../Replication package/';
cd(dir);
addpath(genpath(strcat(dir,"Code/src")))
addpath(genpath(strcat(dir,"Datasets")))

fileID = fopen(strcat(dir,'appl_RUR_data_cleaning.txt'), 'w');

load('di_test');
load('var_names'); 

   n_obs = 388;
   d_var = 131;
  data_m = reshape(data,[n_obs, d_var]);

 d_tilde = data_m(:,varnames=='hlimpw');
 a_tilde = data_m(:,varnames=='accept');   
     bvd = data_m(:,varnames=='bvd');
bvd_a1d1 = (bvd==1);
bvd_a1d0 = (bvd==2);
bvd_a0d1 = (bvd==3);
bvd_a0d0 = (bvd==4);

   white = data_m(:,varnames=='white');
 married = data_m(:,varnames=='married');   
profbavc = data_m(:,varnames=='profbavc');
    male = data_m(:,varnames=='male');
 dispage = data_m(:,varnames=='diapage');

w_1 = [white,married,profbavc,male,dispage];    

  rinc_c = data_m(:,varnames=='rinc_c')/1000;  % income is divided by 1000
   hosps = data_m(:,varnames=='hosps');
 dvisits = data_m(:,varnames=='dvisits');

w_2 = [rinc_c,hosps,dvisits];    
   
  stroke = data_m(:,varnames=='stroke');
 psych_l = data_m(:,varnames=='psych_l');
 arthrit = data_m(:,varnames=='arthrit');
fracture = data_m(:,varnames=='fracture');
backprob = data_m(:,varnames=='backprob');

w_3 = [stroke,psych_l,arthrit,fracture,backprob];

d_walkrm = data_m(:,varnames=='d_walkrm');
   d_sit = data_m(:,varnames=='d_sit');
 d_getup = data_m(:,varnames=='d_getup');
   d_bed = data_m(:,varnames=='d_bed');
d_stairs = data_m(:,varnames=='d_stairs');
dif_eadr = data_m(:,varnames=='dif_eadr');

w_4 = [d_walkrm,d_sit,d_getup,d_bed,d_stairs,dif_eadr];

 wkdpre = data_m(:,varnames=='wkdpre');
hwkdpst = data_m(:,varnames=='hwkdpst');

w_5 = [wkdpre,hwkdpst];

w_matrix = [w_1, w_2, w_3, w_4, w_5];
data_matrix = [a_tilde, d_tilde, w_matrix];

    data_e = rmmissing(data_matrix);
[N, p_all] = size(data_e);
     dataw = zscore(data_e(:,3:end));
     datau = data_e(:,1)-data_e(:,2);
         p = p_all-2;
    accept = data_e(:,1);
    hlimpw = data_e(:,2);

w_var_names = {'white','married','training','male','age',...
               'income','hospital','doctor visits','stroke','psych',...
               'arthrit','fracture','backprob','walking','sitting',...
               'getup','bed','stairs','eat dress','worked prev','hours worked'};             
           
sum_stat = [mean(data_e).' std(data_e).' min(data_e).'  max(data_e).'];

tb = fc_lib_latex_tab_num(sum_stat,1,2,'r',1,0);


% save results
save(strcat(dir,'Datasets/BBCCR.mat'),'dataw','datau','accept','hlimpw','N','p','w_var_names') 

fprintf(fileID,...
    'The sample size is %d and the number of covariates is %d.\n',N,p);

fprintf(fileID,...
    'Table 1. Cross Tabulation: {11,12,21,22} = {a0d0, a0d1, a1d0, a1d1} \n');

fprintf(fileID,...
   '%10d %10d\n%10d %10d\n',transpose(crosstab(data_e(:,1),data_e(:,2))));

fprintf(fileID,...
    'Summary Statistics: Mean\n');

fprintf(fileID,...
    '%8.3f\n',transpose(mean(data_e)));

fprintf(fileID,...
    'Summary Statistics: Median\n');

fprintf(fileID,...
    '%8.3f\n',transpose(median(data_e)));

fprintf(fileID,...
    'Summary Statistics: STD\n');

fprintf(fileID,...
    '%8.3f\n',transpose(std(data_e)));

fprintf(fileID,...
    'Summary Statistics: Min\n');

fprintf(fileID,...
    '%8.3f\n',transpose(min(data_e)));

fprintf(fileID,...
    'Summary Statistics: Max\n');

fprintf(fileID,...
    '%8.3f\n',transpose(max(data_e)));

fprintf(fileID,...
    '%8.3f\n',transpose(max(data_e)));

fclose(fileID);